
# ------------------------------------------------------------------------------------------
# # Julia中的线性代数
# Based on work by Andreas Noack Jensen
#
# ## 基础线性代数操作
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 首先定义一个随机矩阵
# ------------------------------------------------------------------------------------------

A = rand(1:4,3,3)

# ------------------------------------------------------------------------------------------
# 定义一个元素全为1的向量
# ------------------------------------------------------------------------------------------

x = fill(1.0, (3))

# ------------------------------------------------------------------------------------------
# 注意$A$的类型为Array{Int64,2}，而$x$的类型为Array{Float64,1}。Julia定义向量Vector{Type}的别名为Array{Type,1}，矩阵
# Matrix{Type}的别名为Array{Type,2}。
#
# 许多基础操作和其他语言一样
# #### 乘法
# ------------------------------------------------------------------------------------------

b = A*x

# ------------------------------------------------------------------------------------------
# #### 转置
# 就像在其他语言中`A'`表示对`A`进行共轭转置
# ------------------------------------------------------------------------------------------

A'

# ------------------------------------------------------------------------------------------
# 我们可以通过transpose获得转置矩阵
# ------------------------------------------------------------------------------------------

transpose(A)

# ------------------------------------------------------------------------------------------
# #### 转置的乘法
# Julia中某些情况下可以省略`*`号
# ------------------------------------------------------------------------------------------

A'A

# ------------------------------------------------------------------------------------------
# #### 解线性方程组
# 用方阵$A$表示的线性方程组$Ax=b$用左除运算符（函数）`\`求解
# ------------------------------------------------------------------------------------------

A\b

# ------------------------------------------------------------------------------------------
# ## 特殊的矩阵结构
#
# 矩阵结构在线性代数中非常重要。接触一下大一些的线型系统就可以看到矩阵结构有*多*重要了。用线性代数标准包LinearAlgebra可以获得结构化的矩阵structured
# matrices：
# ------------------------------------------------------------------------------------------

using LinearAlgebra

n = 1000
A = randn(n,n);

# ------------------------------------------------------------------------------------------
# Julia可以推断特殊矩阵结构，比如判断对称矩阵
# ------------------------------------------------------------------------------------------

Asym = A + A'
issymmetric(Asym)

# ------------------------------------------------------------------------------------------
# 但有时候浮点错误会比较麻烦
# ------------------------------------------------------------------------------------------

Asym_noisy = copy(Asym)
Asym_noisy[1,2] += 5eps()

issymmetric(Asym_noisy)

# ------------------------------------------------------------------------------------------
# 幸运的是我们可以通过如`Diagonal`，`Triangular`，`Symmetric`，`Hermitian`，`Tridiagonal`和`SymTridiagonal`这
# 样的函数来明确地创建矩阵
# ------------------------------------------------------------------------------------------

Asym_explicit = Symmetric(Asym_noisy);

# ------------------------------------------------------------------------------------------
# 我们来看看Julia计算`Asym`，`Asym_noisy`和`Asym_explicit`的特征值各要花多少时间
# ------------------------------------------------------------------------------------------

@time eigvals(Asym);

@time eigvals(Asym_noisy);

@time eigvals(Asym_explicit);

# ------------------------------------------------------------------------------------------
#
# 本例中，使用`Symmetric()`处理`Asym_noisy` 后让计算效率提高了约5倍
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### A big problem
# Using the `Tridiagonal` and `SymTridiagonal` types to store tridiagonal matrices makes it
# possible to work with potentially very large tridiagonal problems. The following problem
# would not be possible to solve on a laptop if the matrix had to be stored as a (dense)
# `Matrix` type.
# ------------------------------------------------------------------------------------------

n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)
